#include <stdio.h>
#include <graphics.h>
#include <math.h>


double*solvematrix(double**M,double*P,int N)
{
    double*X=new double[N];
 
    int I,J,L;
    double K,S=0;
    double A[N][N];
    double B[N];
    
    for (J=0; J<N; J++) 
    {
      for (I=0; I<N; I++) 
      {
        A[I][J]=*(((double*)M)+I*200+J);
      }
      B[J]=P[J];
    }

    for (I=0; I<=N-2; I++)
    {
      for (J=I+1; J<=N-1; J++)
      {
        K=A[I][I]/A[I][J];
        for (L=I; L<=N-1; L++) 
        {
          A[L][J]=A[L][J]*K-A[L][I];
        }
        B[J]=B[J]*K-B[I];
      }
    }
    
    for (J=N-1; J>=0; J--) 
    {
      S=0;
      for (I=N-1; I>J; I--) 
      {
        S+=A[I][J]*(*(X+I));
      }
      *(X+J)=(B[J]-S)/A[J][J];
    }
    
    return X;
}

int main() 
{
    int N; 
    N=4; 
            
    double A[200][200];
    double B[200];
    double*X=new double[N];

    double Norm;
    int Step=0;
    
    initwindow(400,300,"Norm(N)");
    setcolor(WHITE); 
    setbkcolor(BLACK);

    setfillstyle(1,YELLOW);
    cleardevice();
    int x0=100; 
    int y0=200;   
    int x1=300; 
    int y1=100;   
    

    moveto(x0,y0);
    lineto(x0,y1);
    moveto(x0,y0);
    lineto(x1,y0);

    moveto(x0-3,y1+5);
    lineto(x0,y1);
    moveto(x0+3,y1+5);
    lineto(x0,y1);
    

    moveto(x1,y0);
    lineto(x1-5,y0-3);
    moveto(x1-5,y0+3);
    lineto(x1,y0);

    double E=0.001;

    int J,I;
    J=0;

    X=solvematrix((double**)A,(double*)B,N);
    
    for (int I=0; I<=N-1; I++)
    {
      printf("x%d = %f; ",I,*(X+I));
    }
    printf("\n");
    
    N=1;
    while(true==true)
    {
      N++;
      double * X1=new double[N];
      printf("N = %d\n",N);
      for (J=0; J<=N-1; J++)
      {
        for (I=0; I<=N-1; I++)
        {
          if(I==J)  
          {
            A[I][J]=1+E;
          }
          else if(I<J)
          {
            A[I][J]=E;
          }
          else if(I>J)
          {
            A[I][J]=1-E;
          }
        }
        if(J==N-1)
        {
          B[J]=1+E;
        }
        else
        {
          B[J]=-1;
        }
      }   
      
      for (J=0; J<=N-1; J++) 
      {
        for (I=0; I<=N-1; I++) 
        {
          printf("%f ",A[I][J]);
        }
          printf("| %f\n",B[J]);
      }   

      X1=solvematrix((double**)A,(double*)B,N);
      for (int I=0; I<=N-1; I++)
      {
        printf("x%d = %.15f; ",I,*(X1+I));
        Norm+=X1[I]*X1[I];
      }
      printf("\n");
        
      Norm=pow(Norm,0.5)/N;
      printf("Norm= %.15f\n",Norm);
      fillellipse(101+Step*5,200-(int)(Norm*100),1,1);
      Step++;
      delete X1;
      if(Step>=16)
      {
        break;   
      }
    }  
    
    getchar();
    closegraph();
	return 0;
}
